setwd("...")
options("modelsummary_format_numeric_latex" = "plain")
# options(modelsummary_factory_default = 'kableExtra')
library(modelsummary)
library(kableExtra)
library(ggplot2)
library(broom)
library(paletteer)
library(tidyr)
library(dplyr)
library(stringr)
source("extra/func.R")
theme_set(my_theme())

### ====================================================================

# Load data and limit to days_control_republicano > 0
data = read.csv("data/dataset.csv")
data = subset(data, dias_control_rep > 0)

# Exclude provinces without variation
p_no_variation = data %>%
  group_by(prov_name) %>%
  summarize(ac_bin = paste(unique(anticlerical_bin), collapse = ";")) %>%
  filter(!grepl(";", ac_bin)) %>%
  pull(prov_name)

data = subset(data, !prov_name %in% p_no_variation)

### ====================================================================

## MODEL SPECS

## DVs
dvm = c("anticlerical_muerte_l", "anticlerical_secular_muerte_l", "anticlerical_regular_muerte_l")
dvv = c("anticlerical_l", "anticlerical_secular_l", "anticlerical_regular_l")
dv = c(dvv, dvm)

## Formulae
# Base model
f0 = "patr33_bin + CNT_bin + UGT_bin"
f0_l = "socios_patr33_l + CNT_1936_l + UGT_1932_l"
f0_p_l = "socios_patr33_p_l + CNT_1936_p_l + UGT_1932_p_l"
# Interactions
fint_CNT = "patr33_bin * CNT_bin + UGT_bin"
fint_UGT = "patr33_bin * UGT_bin + CNT_bin"
f_right = c(f0, fint_CNT, fint_UGT)

# Controls
controls = c("pop1930_l", "izq1936", "dias_control_rep_l",
  "analfab30_total", "clerics_all_l")
controls_f = paste(c(controls, "factor(prov_name)"), collapse = " + ")

# Constants (smalltown sample)
smalltown_pop = 10000

### ====================================================================

## MAPS & CIA for TABLES

coef_recode = c(
  "patr33_bin" = "Employers' association",
  'CNT_bin' = 'CNT union',
  'UGT_bin' = 'UGT union',
  'patr33_bin:CNT_bin' = 'Employer $\\times$ CNT union',
  'patr33_bin:UGT_bin' = 'Employer $\\times$ UGT union',
  'pop1930_l' = 'Population (log)',
  'izq1936' = 'Leftist support 1936',
  'compet1936' = 'Electoral competition 1936',
  'dias_control_rep_l' = 'Days of Republican control (log)',
  'analfab30_total' = 'Illiteracy 1930',
  'clerics_all_l' = 'Number of clerics (log)',
  "socios_patr33_l" = "Employers' association (logged)",
  'CNT_1936_l' = 'CNT union (logged)',
  'UGT_1932_l' = 'UGT union (logged)',
  "socios_patr33_p_l" = "Employers' association (pop-weighted, logged)",
  'CNT_1936_p_l' = 'CNT union (pop-weighted, logged)',
  'UGT_1932_p_l' = 'UGT union (pop-weighted, logged)',
  'socios_patr33_l:CNT_1936_l' = 'Employer $\\times$ CNT union',
  'socios_patr33_l:UGT_1936_l' = 'Employer $\\times$ UGT union',
  'socios_patr33_p_l:CNT_1936_p_l' = 'Employer $\\times$ CNT union',
  'socios_patr33_p_l:UGT_1936_p_l' = 'Employer $\\times$ UGT union',
  'compet_UGT_CNT_noNA' = 'Competition UGT vs CNT',
  'compet_UGTCNT_patr33_noNA' = "Competition trade unions vs employers' assoc"
)

gs = list(
  list("raw" = "nobs", "clean" = "$n$", "fmt" = 0),
  list("raw" = "r.squared", "clean" = "$R^2$", "fmt" = 2),
  list("raw" = "adj.r.squared", "clean" = "Adj. $R^2$", "fmt" = 2))

n = "+ p $<$ 0.1, * p $<$ 0.05, ** p $<$ 0.01, *** p $<$ 0.001. Province FE not shown, clustered SE (province). Only municipalities that were at least one day under Republican control."

### ====================================================================
## FORMULA FOR ROBUSTNESS MODELS (w/ interactions)

rob_f_int = function(dv, iv, ctr){

  if(length(iv) > 1){stop("iv must be character vector of the 3 vars")}
  ivs = unlist(str_split(f0, " \\+ "))

  f_base = paste(ivs, collapse = " + ")
  f_cnt = paste(ivs[1], "*", ivs[2], "+", ivs[3])
  f_ugt = paste(ivs[1], "*", ivs[3], "+", ivs[2])
  
  f = c(
    as.formula(paste(dv, "~", f_base, "+", ctr)),
    as.formula(paste(dv, "~", f_cnt, "+", ctr)),
    as.formula(paste(dv, "~", f_ugt, "+", ctr))
  )

  return(f)
}

### ====================================================================
## ALL BASIC MODELS IN ONE LIST
# Base (1) + interactions (2), Base + interactions secular (3) + Base + interactions regular (3))

mod_list_f = apply(expand.grid(f_right, dvv), 1,
  function(x) as.formula(paste(x[2], "~", x[1], "+", controls_f)))
mod_list = lapply(mod_list_f, function(x)
  lm(x, data = data))
names(mod_list) = sapply(mod_list_f, function(x) as.character(x[2]))

### ====================================================================
## MODELS WITH BASE DV

# Get models from list, add names
lm_base = modnames(mod_list[which(grepl("anticlerical_l", names(mod_list)))])

# Table
modelsummary(
  models = lm_base,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence (logged)\\label{tab:lm_base}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lm_base))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_base.tex")

## Predicted probabilities (interaction models)

nd_CNT = expand.grid(patr33_bin = 0:1, CNT_bin = 0:1)
nd_CNT = addcovs(nd_CNT)
nd_UGT = expand.grid(patr33_bin = 0:1, UGT_bin = 0:1)
nd_UGT = addcovs(nd_UGT)

pp_int = bind_rows(
  pp_lm_int(newdf = nd_CNT,
    model = select_mod(mod_list, dv = "anticlerical_l", in_f = "\\* CNT_bin")),
  pp_lm_int(newdf = nd_UGT,
    model = select_mod(mod_list, dv = "anticlerical_l", in_f = "\\* UGT_bin"))) %>%
  mutate(
    patr33_lab = recode(patr33_bin,
      "0" = "No local employers' association",
      "1" = "Local employers' association"),
    org_lab = recode(org,
      "0" = "No local union",
      "1" = "Local union"))

# Plot
p = ggplot(pp_int, aes(x = reorder(org_lab, desc(org_lab)), y = y, group = patr33_lab, color = patr33_lab)) +
  geom_line(position = position_dodge(1/5)) +
  geom_point(position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr90, ymax = upr90), width = 0, linewidth = 1.1,
    position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr95, ymax = upr95), width = 0, position = position_dodge(1/5)) +
  facet_wrap(~ mod_label, ncol = 3) +
  labs(x = "", y = "Predicted clergy killings (log)\n") +
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_color_paletteer_d(`"wesanderson::Darjeeling1"`)
ggsave("output/predprob_int.pdf", height = 3, width = 6, device = "pdf")

### ====================================================================
## MODELS FOR SECULAR/REGULAR

# Get models from list, add names
lm_secreg = modnames(mod_list[which(grepl("secular|regular", names(mod_list)))])

# Table
modelsummary(
  models = lm_secreg,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence by clergy type\\label{tab:lm_secreg}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lm_secreg))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
add_header_above(c(" " = 1, "Secular clergy" = 3, "Regular clergy" = 3)) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_secreg.tex")

### ====================================================================
## ROBUSTNESS MODELS: LINEAR, ALL TOGETHER

# Specs and model list
rob = expand.grid(dv = dv, iv = c(f0, f0_l, f0_p_l), sample = c("full", "smalltowns"))
rob$f = paste0(rob$dv, "~", rob$iv, "+", controls_f)
rob_m = vector("list", nrow(rob))
for(i in 1:nrow(rob)){
  f = as.formula(rob$f[i])
  if(rob$sample[i] == "full"){
    rob_m[[i]] = lm(f, data = data)
  } else if (rob$sample[i] == "smalltowns"){
    rob_m[[i]] = lm(f, data = subset(data, pop1930 < smalltown_pop))
  } else {print(rob$sample[i]);stop("sample?")}
}

### ====================================================================
## ROBUSTNESS: PLOT SUMMARIZING ALL

# Get coefficients
rob$patr33_est = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("patr33", term)) %>% pull(estimate))
rob$patr33_se = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("patr33", term)) %>% pull(std.error))
rob$CNT_est = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("CNT", term)) %>% pull(estimate))
rob$CNT_se = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("CNT", term)) %>% pull(std.error))
rob$UGT_est = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("UGT", term)) %>% pull(estimate))
rob$UGT_se = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("UGT", term)) %>% pull(std.error))
rob$patr33_pval = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("patr33", term)) %>% pull(p.value))
rob$CNT_pval = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("CNT", term)) %>% pull(p.value))
rob$UGT_pval = sapply(rob_m, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("UGT", term)) %>% pull(p.value))

# DF to plot
rob_p = rob %>%
  pivot_longer(
    cols = c(ends_with("_est"), ends_with("_se"), ends_with("_pval")),
    names_to = c("variable", "stat"),
    names_pattern = "^(.*)_(est|se|pval)") %>%
  pivot_wider(names_from = "stat", values_from = "value") %>%
  mutate(
    upr95 = est + qnorm(0.975) * se,
    upr90 = est + qnorm(0.95) * se,
    lwr95 = est - qnorm(0.975) * se,
    lwr90 = est - qnorm(0.95) * se) %>%
  mutate(dv_type = ifelse(grepl("regular", dv), "Regular", "All clergy")) %>%
  mutate(dv_type = ifelse(grepl("secular", dv), "Secular", dv_type)) %>%
  mutate(dv_code = ifelse(grepl("muerte", dv),
    "DV by death place", "DV by residence")) %>%
  mutate(iv_type = ifelse(grepl("_bin", iv), "binary IVs", "continuous IVs (log)")) %>%
  mutate(iv_type = ifelse(grepl("_p_l +", iv),
    "continuous IVs (pop-w and log)", iv_type)) %>%
  mutate(sample_lab = recode(sample,
    "full" = "", "smalltowns" = "excl. > 10000 hab")) %>%
  mutate(model_lab = paste(dv_code, iv_type, sample_lab, sep = "; ")) %>%
  mutate(model_lab = gsub("; $", "", model_lab)) %>%
  mutate(var_lab = recode(variable, "patr33" = "Employers' assoc",
    "CNT" = "CNT unions", "UGT" = "UGT unions", "sind33" = "Trade unions"))

# Plot
p = ggplot(rob_p, aes(x = stringr::str_wrap(model_lab, 50), y = est, color = pval < .05)) +
  geom_pointrange(aes(ymin = lwr90, ymax = upr90), fatten = .5, linewidth = 1.1) +
  geom_pointrange(aes(ymin = lwr95, ymax = upr95), fatten = .75) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  facet_grid(dv_type ~ var_lab) +#, scales = "free") +
  labs(x = "", y = "\nEstimate + 90/95 CIs") +
  coord_flip() +
  theme(
    axis.text.y = element_text(size = 7),
    # strip.text.x = element_text(size = 10, face = "plain"),
    # strip.text.y = element_text(size = 12, face = "plain", angle = 0, hjust = 0),
    # strip.background.y = element_rect(fill = "white", color = NA),
    legend.position = "none") +
  scale_color_manual(values = c("gray", "black"))
ggsave("output/robustness_main.pdf", width = 7, height = 9, device = "pdf")

### ====================================================================
## ROBUSTNESS TABLES: ORG VARIABLES IN CONTINUOUS FORM (builds from previous models)

m_org_l_f = rob_f_int(dv = "anticlerical_l", iv = f0_l, ctr = controls_f)
m_org_l = modnames(lapply(m_org_l_f, function(x) lm(x, data = data)))

m_org_p_l_f = rob_f_int(dv = "anticlerical_l", iv = f0_p_l, ctr = controls_f)
m_org_p_l = modnames(lapply(m_org_l_f, function(x) lm(x, data = data)))

modelsummary(
  models = m_org_l,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence (logged), organization variables in continuous form (logged members)\\label{tab:lm_base_org_l}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_org_l))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_base_org_l.tex")

modelsummary(
  models = m_org_p_l,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence (logged), organization variables in continuous form (logged members)\\label{tab:lm_base_org_p_l}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_org_p_l))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_base_org_p_l.tex")

### ====================================================================
## ROBUSTNESS TABLES: BY PLACE OF DEATH

m_muerte_f = rob_f_int(dv = "anticlerical_muerte_l", iv = f0, ctr = controls_f)
m_muerte = modnames(lapply(m_muerte_f, function(x) lm(x, data = data)))

modelsummary(
  models = m_muerte,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence (logged, by place of death)\\label{tab:lm_base_muerte}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_muerte))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_base_muerte.tex")

### ====================================================================
## ROBUSTNESS TABLES: ONLY SMALL TOWNS

m_smalltowns_f = rob_f_int(dv = "anticlerical_l", iv = f0, ctr = controls_f)
m_smalltowns = modnames(lapply(m_smalltowns_f, function(x)
    lm(x, data = subset(data, pop1930 < smalltown_pop))))

modelsummary(
  models = m_smalltowns,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence (logged), excluding cities (over 10,000 hab.)\\label{tab:lm_base_smalltowns}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_smalltowns))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_base_smalltowns.tex")

### ====================================================================
## ROBUSTNESS TABLES: CONTROLLING FOR COMPETITION (ELECTORAL)

m_elec_compet_f = rob_f_int(dv = "anticlerical_l", iv = f0,
  ctr = gsub("izq1936", "compet1936", controls_f))
m_elec_compet = modnames(lapply(m_elec_compet_f, function(x) lm(x, data = data)))

modelsummary(
  models = m_elec_compet,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence\\label{tab:lm_elec_compet}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_elec_compet))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_elec_compet.tex")

### ====================================================================
## ROBUSTNESS TABLES: CONTROLLING FOR COMPETITION (UGT vs CNT)

m_union_compet_vars = c("compet_UGT_CNT_noNA", "compet_UGTCNT_patr33_noNA")
m_union_compet_f0 = paste("anticlerical_l ~ ", f0, "+ COMPET +", controls_f)
m_union_compet_f = lapply(m_union_compet_vars, function(x)
  as.formula(gsub("COMPET", x, m_union_compet_f0)))

m_union_compet = modnames(lapply(m_union_compet_f, function(x) lm(x, data = data)))

modelsummary(
  models = m_union_compet,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on anticlerical violence\\label{tab:lm_union_compet}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(m_union_compet))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_union_compet.tex")

### ====================================================================
## EXTENSIONS: COMPARE ANTICLERICAL VIOLENCE WITH ALL REPUBLICAN VIOLENCE

data$repub_violence_l_v2 = data$repub_violence - data$anticlerical
data$repub_violence_l_v2 = ifelse(data$repub_violence_l_v2 < 0, 0, data$repub_violence_l_v2)
data$repub_violence_l_v2 = log(data$repub_violence_l_v2 + 1)

# Antic vilence but only with observations where we have data on Repub violence
lm_rep1 = lm(as.formula(paste0("anticlerical_l", "~", f0, "+", controls_f)),
  data = subset(data, !is.na(repub_violence_l_v2)))
# Repub violence (excl anticle violence)
lm_rep2 = lm(as.formula(paste0("repub_violence_l_v2", "~", f0, "+", controls_f)),
  data = data)

# Also including base model for comparison
lm_repub = list(lm_base[[1]], lm_rep1, lm_rep2)
names(lm_repub) = paste0("(", 1:length(lm_repub), ")")

note = paste0(n, " Control variables not shown: ",
    tolower(paste(recode(controls, !!!coef_recode), collapse = ", ")),
  ". Model 2 is the same as Model 1 but restricts the sample to provinces where there is data available on Republican violence (",
  paste(capitalize(unique(data$prov_name[!is.na(data$repub_violence)])), collapse = ", "),
  ").")

modelsummary(
  models = lm_repub,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_rename = coef_recode,
  coef_omit = paste(c(controls, "prov_name"), collapse = "|"),
  gof_map = gs,
  title = "Linear models on Republican violence (logged)\\label{tab:lm_repub}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lm_repub))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
add_header_above(c(" " = 1, "Clerics" = 2, "All victims" = 1)) %>%
add_header_above(c(" " = 1, "Full sample" = 1, "Reduced sample" = 2)) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = note, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_repub.tex")

### ====================================================================
## EXTENSIONS: ALTERNATIVES TO PATRONALES VARIABLE

alt_iv = c(
  "patr33_agra_bin + UGT_bin + CNT_bin +",
  "patr33_no_agra_bin + UGT_bin + CNT_bin +",
  "patr33_agra_bin * UGT_bin + CNT_bin +",
  "patr33_no_agra_bin * UGT_bin + CNT_bin +",
  "patr33_agra_bin * CNT_bin + UGT_bin +",
  "patr33_no_agra_bin * CNT_bin + UGT_bin +"
)

lm_alt_iv = lapply(alt_iv, function(x) {
  lm(as.formula(paste("anticlerical_l ~", x, controls_f)), data = data)
  })
names(lm_alt_iv) = paste0("(", 1:length(lm_alt_iv), ")")

note2 = paste0(n, " Control variables not shown: ",
    tolower(paste(recode(controls, !!!coef_recode), collapse = ", ")), ".")

coef_recode_alt = c(
  "patr33_agra_bin" = "Agrarian employers' assoc",
  "patr33_no_agra_bin" = "Non-agrarian employers' assoc",
  'CNT_bin' = 'CNT union',
  'UGT_bin' = 'UGT union',
  'patr33_agra_bin:CNT_bin' = 'Employer (agrarian) $\\times$ CNT union',
  'patr33_agra_bin:UGT_bin' = 'Employer (agrarian) $\\times$ UGT union',
  'patr33_no_agra_bin:CNT_bin' = 'Employer (not agr.) $\\times$ CNT union',
  'patr33_no_agra_bin:UGT_bin' = 'Employer (not agr.) $\\times$ UGT union')

modelsummary(
  models = lm_alt_iv,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode_alt,
  coef_omit = paste(c(controls, "prov_name"), collapse = "|"),
  gof_map = gs,
  title = "Linear models on anticlerical violence, distinguishing agrarian and non-agrarian employers' associations\\label{tab:lm_alt_iv}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(lm_alt_iv))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = note2, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_agrarian.tex")

# PLOT: compare coefficient of patronales, agrarian, and non-agrarian

ml_patr_agra = list(lm_base[[1]],
  lm_alt_iv[[which(alt_iv == c("patr33_agra_bin + UGT_bin + CNT_bin +"))]],
  lm_alt_iv[[which(alt_iv == c("patr33_no_agra_bin + UGT_bin + CNT_bin +"))]])

coef_patr_agra = lapply(ml_patr_agra, function(x)
  tidy(x, conf.int = TRUE) %>% filter(grepl("patr33", term))) %>%
  bind_rows %>%
  select(term, est=estimate, se=std.error) %>%
  mutate(
    upr95 = est + qnorm(0.975) * se,
    upr90 = est + qnorm(0.95) * se,
    lwr95 = est - qnorm(0.975) * se,
    lwr90 = est - qnorm(0.95) * se) %>%
  mutate(var = recode(term,
    "patr33_bin" = "(Any employers' association)",
    "patr33_agra_bin" = "Agrarian employers' association",
    "patr33_no_agra_bin" = "Non-agrarian employers' association")) %>%
  mutate(all = factor(ifelse(grepl("Any", var), 1, 0)))

# Plot
p = ggplot(coef_patr_agra, aes(x = var, y = est, color = all)) +
  geom_point(position = position_dodge(1/5)) +
    geom_hline(yintercept = 0, linetype = "dotted") +
  geom_errorbar(aes(ymin = lwr90, ymax = upr90), width = 0, linewidth = 1.1,
    position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr95, ymax = upr95), width = 0, position = position_dodge(1/5)) +
  labs(x = "", y = "\nEstimate + 90/95% CIs") +
  theme(legend.title = element_blank(), legend.position = "none") +
  scale_color_manual(values = c("black", "gray")) +
  coord_flip()
ggsave("output/coefplot_patronales_agrarias.pdf",
  height = 3, width = 6, device = "pdf")

### ====================================================================
## EXTENSIONS: TESTING THE ANTICLERICAL TRADITION ISSUE

coef_recode_anticletrad = c(
  "anticle_1820_mc" = "AC violence in Catalonia, 1820s",
  "anticle_zgz_1830s" = "AC violence in Zaragoza, 1830-40s",
  "anticle_2rep_lv" = "AC violence, 1931-1936")

# ----------------------
# Aquillué - Armas y votos: Politización y conflictividad politica en España, 1833-1843
# episodios de 

anticle_zgz_1830s = c(
  50027, # Ambel
  50034, # Ariza
  50037, # Atea
  50074, # Caspe
  50081, # Cetina
  50110, # (Inogés, que se integró en El Fresno)
  50137, # Leciñena
  50143, # Longares
  50200, # Paniza
  50297, # Peñaflor (integrado en Zaragoza en 1890s)
  50209, # Pinseque
  50252, # Tauste
  50293, # Villaroya (imagino será V. de la Sierra)
  50297) # Zaragoza

data$anticle_zgz_1830s = ifelse(data$COD_INE %in% anticle_zgz_1830s, 1, 0)

# Base formula
base = paste("anticlerical_l ~",
  paste(f0, paste(controls, collapse = " + "), sep = " + "))

# Models, with/without variables
# Anticlerical violence Catalunya in 1822-23, Muns y Casteller (1888)
lm_mc1 = lm(as.formula(paste(base, "+ factor(prov_name)")),
  data = subset(data, ccaa == "catalunya"))
lm_mc2 = lm(as.formula(paste(base, "+ factor(prov_name) + anticle_1820_mc")),
  data = subset(data, ccaa == "catalunya"))
# Violencia AC, 1834-1843 en Zaragoza (Aquillue)
lm_zgz1 = lm(as.formula(base),
  data = subset(data, prov_name == "zaragoza"))
lm_zgz2 = lm(as.formula(paste(base, "+ anticle_zgz_1830s")),
  data = subset(data, prov_name == "zaragoza"))
# Lopez-Villaverde, motines anticlericales Segunda Republica
lm_lv1 = lm(as.formula(paste(base, "+ factor(prov_name)")),
  data = data)
lm_lv2 = lm(as.formula(paste(base, "+ factor(prov_name) + anticle_2rep_lv")),
  data = data)

# Model list
ml_anticletrad = list(lm_mc1, lm_mc2, lm_zgz1, lm_zgz2, lm_lv1, lm_lv2)

# Note for table
n_trad = paste(gsub(", clustered SE \\(province\\)", "", note2),
  "Model 2 controls for anticlerical violence in Catalonia in 1822-23, following Muns i Casteller (1888). Model 4 controls for the incidence of anticlerical violence in Zaragoza province between 1833 and 1843 following Aquillué (2020). Model 6 controls for the incidence of anticlerical mutinies during the Second Republic, 1931-1936, following López Villaverde (2019).")

# Table
modelsummary(
  models = ml_anticletrad,
  output = "latex",
  # vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = c(coef_recode, coef_recode_anticletrad),
  coef_omit = paste(c(controls, "prov_name"), collapse = "|"),
  gof_map = gs,
  title = "Linear models on anticlerical violence controlling for historical anticlericalism\\label{tab:lm_trad_anticle}",
  add_rows = as.data.frame(rbind(c("Province FE", rep("Yes", length(ml_anticletrad))))),
  threeparttable = TRUE,
  escape = FALSE) %>%
add_header_above(c(" " = 1, "Catalonia" = 2, "Zaragoza" = 2, "Full sample" = 2)) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n_trad, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_trad_anticle.tex")

### ====================================================================
## EXTENSIONS: ONLY NUNS

nuns_provs_t = table(data$prov_name, data$anticlerical_mujeres_bin)
nuns_provs = rownames(nuns_provs_t)[nuns_provs_t[,2] > 0]

lm_nuns1 = lm(as.formula(paste("anticlerical_mujeres_l ~",
  paste(f0, paste(controls_f, collapse = " + "), sep = " + "))),
  data = data)
lm_nuns2 = glm(as.formula(paste("anticlerical_mujeres_bin ~",
  paste(f0, paste(controls_f, collapse = " + "), sep = " + "))),
  data = subset(data, prov_name %in% nuns_provs), family = "binomial")
lm_nuns3 = glm(as.formula(paste("anticlerical_mujeres_bin ~",
  paste(f0, paste(controls, collapse = " + "), sep = " + "))),
  data = data, family = "binomial")

ml_nuns = list(lm_nuns1, lm_nuns2, lm_nuns3)

n_nuns = paste(n, "Model 1 includes the logged number of female clerics killed as dependent variable, and Models 2 and 3 use a binary indicator as dependent variable. According to our data, we only register killings of female clerics in 48 municipalities (slightly over 1 per cent). Due to perfect separation problems, Model 2 includes province FE but limits the sample to provinces with variation, while Model 3 includes the full sample without province FE.")


gs_nuns = list(
  list("raw" = "nobs", "clean" = "$n$", "fmt" = 0),
  list("raw" = "r.squared", "clean" = "$R^2$", "fmt" = 2),
  list("raw" = "adj.r.squared", "clean" = "Adj. $R^2$", "fmt" = 2),
  list("raw" = "aic", "clean" = "AIC", "fmt" = 2))

# Table
modelsummary(
  models = ml_nuns,
  output = "latex",
  vcov = ~prov_name,
  estimate = "{estimate}{stars}",
  coef_map = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs_nuns,
  title = "Linear models on anticlerical violence against female clerics (nuns)\\label{tab:lm_nuns}",
  add_rows = as.data.frame(rbind(c("Province FE", c("Yes", "Yes", "No")))),
  threeparttable = TRUE,
  escape = FALSE) %>%
add_header_above(c(" " = 1, "LM (log)" = 1, "Logit (binary)" = 2)) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n_nuns, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_nuns.tex")
